rm(list = ls(all = TRUE))
library(foreign)
library(lme4)
library(Amelia)
library(plyr)
library(texreg)

GG <- read.dta("full_parliaments.dta")

## Functions
aggrMI <- function(mod){
    out <- list(NULL)
    beta <- matrix(NA, nrow = 5, ncol = 3)
    sigm <- matrix(NA, nrow = 5, ncol = 3)
    for(i in 1:5){
        beta[i, ] <- getME(mod[[i]], name = "beta")
        sigm[i, ] <- as.matrix(diag(vcov(mod[[i]])))
    }
    
    betabar <- round(colMeans(beta), digits = 3)
    outses <- rep(NA, 3)
    m <- 5
    for(i in 1:3){## See p. 53 in King et. al. 2001
        a <- sum((beta[, i] - betabar[i])^2) / (m - 1)
        b <- sum(sigm[, i]) / m
        c <- b + (a * (1 + (1 / m)))
        outses[i] <- sqrt(c)
    }
    out[[1]] <- betabar
    out[[2]] <- outses
    return(out)
}

## Table 4
GGH <- GG[GG$deposed=="P", ]
GGH$name_country <- paste(GGH$name, GGH$country, GGH$duration)
GGH$presence <- as.numeric(as.character(GGH$presence))
GGH$duration <- as.numeric(as.character(GGH$duration))
GGH$age <- as.numeric(as.character(GGH$age))
ggh <- ddply(GGH, c("name_country"), summarise, d = sum(presence), n = length(presence),
               duration = unique(duration), age = unique(age), country = unique(country),
             house = unique(house), century = unique(century), name = unique(name))

ggh$house <- as.numeric(ggh$house)
ggh$century <- as.numeric(ggh$century)
ggh$country <- as.numeric(ggh$country)
ggh$lnage <- log(ggh$age + 1)
ggh$lnage <- ifelse(is.na(ggh$lnage), NA, ggh$lnage)
ggh$lnduration <- log(ggh$duration + 1)


binom1 <- glm(I(d / n) ~ lnduration, weights = n, family = binomial(link=logit),data = ggh)
binom2 <- glmer(I(d / n) ~ lnduration + (1 | country), weights = n, family = binomial(link=logit), data = ggh)
binom3 <- glmer(I(d / n) ~ lnduration + (1 | country) + (1 | century), weights = n, family = binomial(link=logit), data = ggh)
binom4 <- glmer(I(d / n) ~ lnduration + (1 | country) + (1 | century) + (1 | house), weights = n, family = binomial(link=logit), data = ggh)
set.seed(444)
a.out <- amelia(ggh[, c("lnage", "country","house",  "century", "lnduration","name","n","d")], m = 5, ts = c("century"), cs = c("name"))
fA <- formula(I(d / n) ~ lnduration + lnage + (1 | country) + (1 | century ) + (1 | house))
binom5 <- lapply(a.out$imputations, function(x) glmer(formula = fA, weights = n, family = binomial(link=logit), data = x))
binom5aggr <- aggrMI(binom5)

beta <- list(coef(binom1), getME(binom2, name = "beta"), getME(binom3, name = "beta"), getME(binom4, name = "beta"), binom5aggr[[1]])
ses <- list(as.numeric(sqrt(diag(vcov(binom1)))), sqrt(diag(vcov(binom2))), sqrt(diag(vcov(binom3))), sqrt(diag(vcov(binom4))), binom5aggr[[2]])
pval <- list(NULL)
for(i in 1:length(beta)){
    temp.beta <- beta[[i]]
    temp.ses <- ses[[i]]
    temp.pval <- 2 * (1 - pnorm(abs(temp.beta / temp.ses)))
    pval[[i]] <- as.numeric(temp.pval)
}

texreg(list(binom1, binom2, binom3, binom4, binom5[[1]]), override.coef = beta, override.se = ses, override.pval = pval,
       stars = c(0.01, 0.05, 0.1), digits = 4, file = "table4.txt")

save(binom1, file = "parliament_binom.RData")
